///@file Tensor2.cpp
///@brief Implementation of the class representing a second order tensor
///@author Arnaud Duval
///@version 0.0
///@date 2009/10/29


#include <assert.h>
#include <cmath>

#include "Tensor2.h"

/*///@brief Desctuctor
Tensor2::~Tensor2()
{
	delete[] arows;
}*/

///@brief Return dimension
unsigned int Tensor2::Dim()
{
	assert(nbRows == nbColumns);
	return nbRows;
}

///@brief Return deviatoric part of tensor
Tensor2 Tensor2::Dev()
{
	assert(nbColumns == nbRows);
	assert(nbColumns != 0);
	Tensor2 res(nbColumns);
	double trace = 0.;
	for(unsigned int i = 1 ; i <= nbColumns ; i++)
		trace += (*this)(i,i);
	
	for(unsigned int i = 1 ; i <= nbColumns ; i++)
		for(unsigned int j = 1 ; j <= nbColumns ; j++)
		{
			if(i == j)
				res(i,j) = (*this)(i,j) - trace/nbColumns;
			else
				res(i,j) = (*this)(i,j);
		}
	
	return res;
}

///@brief Return a Kronecker tensor
///@param dim Dimension of the tensor ti return
Tensor2 Tensor2::Delta(unsigned int dim)
{
    Tensor2 res(dim);
    for(unsigned int i = 1 ; i <= dim ; i++)
        for(unsigned int j = 1 ; j <= dim ; j++)
            if(j == i)
                res(i,j) = 1.;
            else
                res(i,j) = 0.;
    return res;
}

///@brief Return Von Mises equivalent for stresses
double Tensor2::VMEquivStress()
{
    Tensor2 d = Dev();
    return sqrt((3./2.)*(d||d));
}

///@brief Overload of + operator
///@param t Tensor2 right hand operand
Tensor2 Tensor2::operator+ (const Tensor2& t)
{
    assert(nbColumns == t.nbColumns);
    Tensor2 res(nbColumns);
    for(unsigned int i = 1 ; i <= nbColumns ; i++)
        for(unsigned int j = 1 ; j <= nbColumns ; j++)
            res.arows[i-1][j] = arows[i-1][j] + t.arows[i-1][j];
    return res;
}
        
///@brief Overload of - operator
///@param t Tensor2 right hand operand
Tensor2 Tensor2::operator- (const Tensor2& t)
{
    assert(nbColumns == t.nbColumns);
    Tensor2 res(nbColumns);
    for(unsigned int i = 1 ; i <= nbColumns ; i++)
        for(unsigned int j = 1 ; j <= nbColumns ; j++)
            res.arows[i-1][j] = arows[i-1][j] - t.arows[i-1][j];
    return res;
}

///@brief Overload of / operator : division by a scalar
///@param s Scalar right hand operand
Tensor2 Tensor2::operator/ (const double& s)
{
    
    Tensor2 res(nbColumns);
    if(s != 0.)
        for(unsigned int i = 1 ; i <= nbColumns ; i++)
            for(unsigned int j = 1 ; j <= nbColumns ; j++)
                res.arows[i-1][j] = arows[i-1][j] / s;
    else if(VMEquivStress() == 0.)
        res.Nullify();
    else
        std::cerr << "Error : division by zero in Tensor2::operator/\n";
    return res;
}

///@brief Overload of / operator : multiplication by a scalar
///@param s Scalar right hand operand
Tensor2 Tensor2::operator* (const double& s)
{
    Tensor2 res(nbColumns);
    for(unsigned int i = 1 ; i <= nbColumns ; i++)
        for(unsigned int j = 1 ; j <= nbColumns ; j++)
            res.arows[i-1][j] = arows[i-1][j] * s;
    return res;
}

///@brief Overload of * operator : product scalar * tensor
Tensor2 operator* (const double& s, const Tensor2& tens)
{
    Tensor2 res(tens.nbColumns);
    for(unsigned int i = 1 ; i <= tens.nbColumns ; i++)
        for(unsigned int j = 1 ; j <= tens.nbColumns ; j++)
            res.arows[i-1][j] = tens.arows[i-1][j] * s;
    return res;
}
        
///@brief Overload of | operator : contracted product A_ij = B_ik C_kj
///@param t1 left hand operand
///@param t2 right hand operand
Tensor2 operator | (const Tensor2& t1, const Tensor2& t2)
{
  assert (t1.nbColumns == t2.nbColumns);
  assert (t1.nbRows == t2.nbRows);
  assert (t1.nbColumns == t1.nbRows);
  Tensor2 res(t1.nbRows);
  for(unsigned int i = 1 ; i <= t1.nbRows ; i++)
  {
    for(unsigned int j = 1 ; j <= t1.nbRows ; j++)
    {
      res.arows[i-1][j] = 0.;
      for(unsigned int k = 1 ; k <= t1.nbRows ; k++)
      {
	res.arows[i-1][j] += t1.arows[i-1][k] * t2.arows[k-1][j];
      }
    }
  }

  return res;
}


///@brief Overload of || operator : double contracted product A = B_ij C_ji
///@param t1 left hand operand
///@param t2 right hand operand
double operator || (const Tensor2& t1, const Tensor2& t2)
{
  assert (t1.nbColumns == t2.nbColumns);
  assert (t1.nbRows == t2.nbRows);
  assert (t1.nbColumns == t1.nbRows);
  double res = 0.;
  for(unsigned int i = 1 ; i <= t1.nbRows ; i++)
      for(unsigned int j = 1 ; j <= t1.nbRows ; j++)
          res += t1.arows[i-1][j] * t2.arows[j-1][i];
  return res;
}

